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Abstract 

Recurrence plot based methods are highly efficient and widely accepted tools for 
the investigation of time series or one-dimensional data. We present an extension of 
the recurrence plots and their quantifications in order to study recurrent structures 
in higher-dimensional spatial data. The capability of this extension is illustrated on 
prototypical 2D models. Next, the tested and proved approach is applied to assess 
the bone structure from CT images of human proximal tibia. We find that the 
spatial structures in trabecular bone become more self-similar during the bone loss 
in osteoporosis. 
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1 Introduction 



Recurrence is a fundamental property of many dynamical systems and, hence, 
of various processes in nature. A system may strongly diverge, but after some 
time it recurs "infinitely many times as close as one wishes to its initial state" 
[1]. The investigation of recurrence reveals typical properties of the system 
and may help to predict its future behaviour. With the study of nonlinear 
chaotic systems several methods for the investigation of recurrences have been 
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developed. The method of recurrence plots (RPs) was introduced by Eck- 
mann et al. [2]. Together with different RP quantification approaches [3,4], 
this method has attracted growing interest for both theory and applications 
[5]. 

Recurrence plot based methods have been succesfully applied to a wide class of 
data from physiology, geology, physics, finances and others. They are especially 
suitable for the investigation of rather short and nonstationary data. This 
approach works with time series or phase space reconstructions (trajectories), 
i. e. with data which are at least one-dimensional. 

Recurrences are not restricted to one- dimensional time series or phase space 
trajectories. Spatio-temporal processes can also exhibit typical recurrent struc- 
tures. However, RPs as introduced in [2] cannot be directly applied to spatial 
(higher-dimensional) data. One possible way to study the recurrences of spa- 
tial data is to separate the higher- dimensional objects into a large number 
of one- dimensional data series, and to analyse them separately [6]. A more 
promising approach is to extend the one-dimensional approach of the recur- 
rence plots to a higher-dimensional one. 

In the presented work, we focus on the analysis of snapshots of spatio-temporal 
processes, e.g., on static images. An extension of recurrence plots and their 
quantification to higher-dimensional data is suggested. This extension allows 
us to apply this method directly to spatial higher-dimensional data, and, in 
particular, to use it for 2D image analysis. We apply this method to 2D pQCT 
human bone images in order to investigate differences in trabecular bone struc- 
tures at different stages of osteoporosis. 



2 Recurrence Plots 

The initial purpose of recurrence plots was the visualisation of recurrences of 
system's states Xi in a phase space (with dimension m) within a small devia- 
tion e [2]. The RP efficiently visuafises recurrences even for high dimensional 
systems. A recurrence of a state at time i at a different time j is marked within 
a two-dimensional squared matrix with ones and zeros dots (black and white 
points in the plot), where both axes represent time. The RP can be formally 
expressed by the matrix 

Rij^e{e-\\xi-xj\\), XiER"", i,j^l...N, (1) 

where is the number of considered states threshold distance (an 

arbitrary deviation range within a recurrence is defined), || • || denotes a norm 
and ©(•) is the Heaviside function. 
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It shoTild be emphasised that this method is a pairwise comparison of system's 
states at different times along a phase space trajectory, which is - affhough 
lying in an m-dimensional space - a one-dimensional curve. The axes of the RP 
correspond to the time which is given by pursueing a state on the trajectory. 
Diagonal lines in an RP represent epochs of similar dynamical evolution of 
the analysed system. For i = j we get the line of identity (LOI), Rj^j = 1 
which is the main diagonal line in the RP (Fig. 1). 

Instead of using the system's states Xi which are often unknown, RPs can be 
created by only using a single time series or a reconstruction of the phase 
space vectors (e.g. by using time-delay embedding, [7]). Such applications to 
experimental data have expanded the utilisation of RPs from a tool for the 
investigation of deterministic phase space dynamics to a tool for the investi- 
gation of similarity and transitions in data series, without the rather strong 
requirement that the data must be from a deterministic dynamical process. 
The idea of such a similarity plot is not new and can be found in publications 
earlier than [2], e.g. in [8]. This alternative understanding was (unconsciously) 
the base of the ever increasing amount of application of RPs in data analy- 
sis. However, in its present state the RP technique could not be applied on 
higher-dimensional spatial data. 



The initial purpose of RPs was the visual inspection of the behaviour of phase 
space trajectories. The appearance of RPs gives hints about the characteristic 
time evolution of these trajectories [5]. A closer inspection of RPs reveals 
small-scale structures which are single dots, diagonal lines as well as vertical 
and horizontal lines (Fig. 1). 




Fig. 1. Example of a recurrence plot for the logistic map (xj+i = axi{l — x) with 
control parameter a = 3.9767). The RP consists of single dots and line structures. 



A diagonal line 'Ri+kj+k = 1 11=0 (where I is the length of the diagonal line) 
occurs when one segment of the trajectory runs parallel to another one, i. e. the 
trajectory re-visits the same region of the phase space at different time inter- 
vals. The length of this diagonal line is determined by the duration of intervals 
with similar local behaviour of the trajectory segments. We define a line in 
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the RP as a diagonal line of length I, if it fulfils the condition 

i-i 

(1 - Rj_ij_i) (1 - Ri+ij+i) Yi R-i+fej+fc = 1- (2) 

fe=0 



A vertical (horizontal) line Rij+jt = 1 I^Iq (where v is the length of the vertical 
line) marks a time interval in which a system's state does not change in time or 
changes very slowly. It looks like the state is trapped for some time, which is a 
typical behaviour of laminar states [4] . Because RPs are symmetric about the 
LOI by definition (1), each vertical line has a corresponding horizontal line. 
Therefore, only the vertical lines are henceforth considered. Combinations of 
vertical and horizontal lines form rectangular clusters in an RP. We define a 
line as a vertical line of length v, if it fulfils the condition 

v+l 

(1 - Ri,j-l) (1 - RjJ+t)) n ^ij+k = 1- (3) 

A:=0 



These small-scale structures are used for the quantitative analysis of RPs 
(known as recurrence quantification analysis, RQA). Using the distributions 
of the lengths of diagonal lines P{1) or vertical lines P{v), different measures 
of complexity have been introduced (cf. [5] for a comprehensive review of def- 
initions and descriptions of these measures) . Here we generalise the measures 
recurrence rate RR, determinism DET, averaged diagonal line length L, 1am- 
inarity LAM and trapping time TT in order to quantify higher-dimensional 
data. (cf. Tab. 1). 

Several measures need a predefined minimal length Imin or Vmin, respectively, 
for the definition of a diagonal or vertical line. These minimal lengths should 
be as minimal as possible in order to cover as much variation of the lengths 
of these lines. On the other hand, Imin a-nd Vmin should be large enough to 
exclude line-like structures which represent only single, non-recurrent states, 
which may occur if the threshold e is chosen too large or if the data have been 
smoothed too stronlgy before computing the RP. 

RQA was successfully applied for example for the detection of transitions in 
event related EEC potentials [9], the study of interrelations between El Nifio 
and climate in the past [10], the investigation of economic data series [11], 
of nonlinear processes in electronic devices [12] or the study of transitions in 
chemical reactions [13]. For a number of further applications see, e.g., [5] or 
www.recurrence-plot.tk. 
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Table 1 

Generalised recurrence quantification measures for spatial data of dimension d and 
with I, J G Z*^. Note that these measures assess recurrence information in terms of 
length while the original RQA measures quantify it in terms of time. 

RQA measure equation meaning 

percentage of recurrent 
states in the system; 
probability of the recur- 
rence of any state 

percentage of recurrence 
points which form diag- 
onal hyper-surfaces; re- 
lated to the predictability 
of the system 

percentage of recurrence 
points which form ver- 
tical hyper-surfaces; re- 
lated to the laminarity of 
the system 

related to the prediction 
length of the system 

related to the size of the 
area in which the system 
does not change 

3 Extension to higher dimensions 

Now, we propose an extension of RPs to analyse higher dimensional data. 
With this step we leave the RPs as a method for investigating deterministic 
dynamics and focus on its potential in determining similar (recurrent) features 
in spatial data. 

For a d-dimensional (Cartesian) system, we define an n-dimensional recurrence 
plot by 

R,,j =Q{£-\\xt- %|| ) , xt e E™, I, J G (4) 

where lis the c?- dimensional coordinate vector and is the phase-space vector 
at the location given by the coordinate vector z. This means that we decom- 
pose the spatial dimension of and consider each space direction separately, 
e. g. .'?ji.i2....id f*^^ ^1 — Ij • • • , ^ but i2, ■ ■ ■ ,id = const. Such vectors are now 
one-dimensional curves in the m-dimensional space. Each of these vectors is 
pairwisely compared with all others. These individual sub-RPs are the com- 
ponents of the final higher-dimensional RP. The resulting RP has now the 



recurrence rate 



iV 



determinism 
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dimension n = 2 x d and cannot be visualised anymore. However, its quantifi- 
cation is still possible. 

Similarly to the one-dimensional LOI given by Hij = 1 \/i — j, we can find 
diagonally oriented, d-dimensional structures in this n- dimensional recurrence 
plot {n — 2d), called the hyper-surface of identity (HSOI): 

Rtj=i yt = j. (5) 



In the special case of a two-dimensional image composed by scalar values, we 
have 

which is in fact a four- dimensional recurrence plot, and its HSOI (R-n, 22,11,12 = 
1) is a two-dimensional plane. 



4 Quantification of Higiier-Dimensional RPs 



The known RQA is based on the quantification of the line structures in the 
two-dimensional RPs. Thus, the definition of higher-dimensional equivalent 
structures is crucial for a quantification analysis of higher-dimensional RPs. 

Based on the definition of diagonal lines, Eq. (2), we define a diagonal squared 
hyper-surface of size I {I — {I, ... ,1), I E Z'*) as 

(1 - Rr-iV-r) (1 - R.w+r) n (7) 

fel,fe2,..., 
fed=0 

In particular, for the two-dimensional case such a diagonal hyper-surface (HS) 
is thus defined as 

l-l 

(1 — Rii_l,i2-l,jl-l,i2-l) (1 ~ ^il+l,i2+l,jl+l,j2+l) n Ril+fcl,i2+fc2,il+fcl,j2+fc2 = ^ 

fel,fe2=0 

(8) 

The next characteristic structure, the vertical squared HS of size v {v — 
{v, . . . ,v),v e Z'^) is defined as 

(1 - Rr,j--„-) (1 - Rr,j-+^;) U ^rj+k^^- (9) 

fel,fe2,..., 
kd=0 
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Its 2D equivalent is 

v-l 

(1 ~ 1^11,12,^1 — 1,^2-1) (-'- ~ ^il,i2,jl+V,j2+v) Y\. ■^n,i2,jl+fel,j2+fe2 = (10) 

fcl,fc2=0 



Using these definitions, we can construct the frequency distributions P{1) and 
P{v) of the sizes of diagonal and vertical HS in the higher-dimensional RP. 
This way we get generalised RQA measures DEThs, LAMhsi Lhs and TThs 
as defined in Tab. 1, which are now suitable for characterising spatial (e. g. two- 
dimensional) data. 



5 Model Examples 



In order to illustrate the ability of the proposed high-dimensional RP's exten- 
sion, we consider three prototypical model examples from 2D image analysis. 
The first image (A) is produced by uniformly distributed white noise, the 
second one (B) is the result of a two-dimensional auto-correlated process of 
2nd order (2D-AR2; ajj ,,• = Yllk 1=1 '^k,i ^i-kj-i + where ttk^i is the 2D matrix 
of model parameters and { is Gaussian white noise) and the third one (C) 
represents periodical recurrent structures (Fig. 2). All these example images 
have a size of 200 x 200 pixels and are normalised to a mean of zero and a 
standard deviation of one. 

The resulting RPs are four-dimensional matrices of size 200 x 200 x 200 x 200, 
and can hardly be visualised. However, in order to visualise these RPs, we 
can reduce their dimension by one by considering only those part of the RPs, 
where ^2 = J2- The resulting 200 x 200 x 200 cube is a hypersurface of the 
four-dimensional RP along the LOI. For the threshold we use e = 0.2, which 
gives clear representations of the RPs. 





Fig. 2. Two-dimensional prototypical examples: test images representing (A) uni- 
formly distributed white noise, (B) a two-dimensional auto-correlated process 
(2D-AR2) and (C) periodical recurrent structures. 
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ABC 

Fig. 3. Three-dimensional subsections of four-dimensional RPs of the images shown 
in Fig. 2. As known from one-dimensional RQA, (A) random data causes homoge- 
neous RPs consisting of single, dis-connected points, (B) correlations in data cause 
extended connected structures and (C) periodic data induce periodically occuring 
structures in the RPs. 



A 




Fig. 4. Slices of the subsections of the four-dimensional RPs shown in Fig. 3. The 
similarity to known recurrence plots is obvious: (A) noise, (B) auto-correlated data 
and (C) periodic data. 

The features occuring in higher-dimensional RPs provide similar information 
as known from the classic one-dimensional RPs. Separated single points cor- 
respond to strongly fluctuating, uncorrelated data as it is typical for, e.g., 
white noise (Fig. 3A). Auto-correlations in data cause extended structures, 
which can be lines, planes or even cuboids (Fig. 3B). Periodical recurrent pat- 
terns in data imply periodic line and plane structures in the RP (Fig. 3C). 
Two-dimensional slices through such RPs contain similar patterns found by 
common RPs (Fig. 4). 

We compute the proposed RQA measures (Tab. 1) for the histograms of the 
sizes of diagonal and vertical planes (2D HS) in the four-dimensional RPs. For 
all three examples we use for the minimal size of the diagonal and vertical 
HS Imin = 3 pixels and Vmin = 4 pixels. Although the RQA measures depend 
on the value of e, its selection is not crucial for our purpose to discriminate 
the three different types of structures in the test images. The chosen values 
for Imin and Vmin are found to be optimal for discriminating the considered 
images. By choosing smaller values of Imin and Vmin (but larger than one), the 
measures DEThs and LAMhs are closer for the 2D-AR2 and the periodic 
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Table 2 

Recurrence quantification measures for the prototypical examples shown in Fig. 2. 
The measures are explained in Tab. 1. 



Example 


RR 


DEThs 


LAMhs 


Lhs 


TThs 


(A) noise 


0.218 


0.007 


0.006 


3.7 


3.0 


(B) 2D-AR2 


0.221 


0.032 


0.065 


3.1 


3.1 


(C) periodic 


0.219 


0.322 


0.312 


5.8 


5.6 



image. 

Four of five RQA measures clearly discriminate between the three types of 
images (Tab. 2). Only the recurrence rate RR is roughly the same for all test 
objects. This is because all images were normahsed to the same standard de- 
viation. For the random image (A) the determinism DEThs and laminarity 
LAMhs tend to zero, what is expected, because the values in the image heav- 
ily fluctuate even between adjecent pixels. For the 2D-AR2 image (B), DEThs 
and LAMhs are slightly above zero, revealing the correlation between adjecent 
pixels. The last example (C) has, as expected, the highest values in DEThs 
and LAMhs, because same structures occur many times in this image and the 
image is rather smooth. Although the trend in DEThs and LAMhs is simi- 
lar, there is a significant difference between both measures. Whereas LAMhs 
represents the probability that a specific value will not change over spatial vari- 
ation (what results in extended same-coloured areas in the image), DEThs 
measures the probability that similar changes in the image recur. LAMhs is 
twice of DEThs for the 2D-AR2 image, obtaining that there are more areas 
without changes in the image than such with typical, recurrent changes. 



6 Application to pQCT data of proximal tibia 

According to the definition of the World Health Organisation, osteoporosis is 
a disease characterised by bone loss and changes in the structure of the bone. 
In the last years, the focus changed to structural assessment of the trabecular 
bone, because bone densitometry alone is very limited to explain all variation 
in bone strength. Furthermore, the rapid progress in the development of new 
high-resolution computer tomography (CT) scanners facilitates investigations 
of the bone micro-architecture. Different approaches using methods coming 
from nonlinear dynamics have been recently proposed in order to evaluate 
structural changes [14,15,16,17] or even to predict fracture risks or biome- 
chanical properties [18,19,20]. These approaches use, e.g., scaling properties 
of bone micro-structure or symbol-encoding of the bone architecture. 

Using the RP based method, we will focus here on the recurrent structures 
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found in images of trabecular bone of proximal tibiae obtained by periph- 
eral quantitative computer tomography (pQCT). The images were aquired 
from bone specimens with different stages of osteoporosis as assessed by bone 
mineral density (BMD). Being applied to such images, the RP provides infor- 
mation about recurrences of bone and soft tissue. 

The spatial recurrence analysis is applied to high-resolution pQCT axial slices 
of human proximal tibia acquired 17 mm below the tibial plateau, with pixel 
size 200 fim and slice thickness 1 mm (Fig. 5). The images were acquired from 
25 bone specimens with a pQCT scanner XCT-2000 (Stratec GmbH, Ger- 
many). The trabecular bone mineral density of these specimens ranges from 
30 to ISOmg/cm"^. A standardised image pre-processing procedure was applied 
to exclude the cortical shell from the analysis [15,21]. The RQA was computed 
by using the parameters e — 0.04 cm~^, — Vmin — 400 /im. These minimal 
lengths correspond to two pixels and is found to be appropriate for pQCT 
images of such resolution. 

In order to further evaluate the proposed RQA measures, we compare them 
with some recently introduced structural measures of complexity (SMCs) 
[15,21]. The SMCs are based on a symbol-encoding of bone elements in the 
pQCT image. Here we focus on the following SMCs: 

(1) Entropy {Sa)- quantifies the probability distribution of X-ray attenuation 
within the ROI; 

(2) Structure Complexity Index (SCI): assesses the complexity and homo- 
geneity of the structure whole; 

(3) Trabecular Network Index (TNI) evaluates richness, orderliness, and ho- 
mogeneity of the trabecular network. 

The computation of the SMCs is applied on the same trabecular area like the 
RQA measures. 




Fig. 5. Typical axial pQCT slice of human proximal tibia acquired 17 mm below 
the tibial plateau. The trabecular BMD is 65.5mg/cm^. 

The application of the recurrence plot extension to the pQCT images of prox- 
imal tibiae reveals a relationship between the recurrences in the trabecular 



10 



architecture and the osteoporotic stage (Fig. 6 and Tab. 3). RR is largest 
for osteoporotic bone and shows the strongest relationship with the degree of 
osteoporosis: it is clearly anti-correlated with BMD (Spearman's rank order 
correlation coefficient R — —0.94). DEThs and LAMhs are also maximal for 
tibiae with high degree of osteoporosis {R — —0.66 and —0.79; Fig. 7). We 
do not find a strong relation between LhS) TThs and BMD. The comparison 
with the SMCs reveals good relationships between the RQA measures and Sa, 
SCI and TNI (Fig. 8 and Tab. 3). Thus, the RQA measures RR, DEThs and 
LAMhs contain also information about the complexity and homogeneity of 
the trabecular network. 

Thus, the proposed RP approach reveals that during the development of os- 
teoporosis the structures in the corresponding pQCT image become more and 
more recurrent. This is in a good agreement with a decreasing complexity in 
the micro-architecture of bone. It confirms the results of an analysis of pQCT 
images acquired from human proximal tibia and lumbar vertebrae based on 
symbolic dynamics [15,21]. The direct comparison with the structural quan- 
tities (SMCs) shows that the RQA measures provide information about the 
bone architecture. The RQA measures reveal a low rate of change for bone 
of higher BMD, but higher rate of changes for specimens with lower BMD 
(Figs. 6 and 7). This reflects a higher sensitivity of these measures for osteo- 
porotic trabecular bone and emphasises the nonlinear relationship between 
the bone architecture as assessed by the RQA measures and bone mass as 
evaluated by the BMD. As it has been recently shown that the SMCs provide 
a better estimation of the mechanical bone strength than BMD alone [22], the 
proposed RQA measures could further enhance the evaluation to assess the 
fracture risk of bone. 



0.5 




50 100 150 

BMD 



Fig. 6. Recurrence rate RR obtained from four-dimensional RPs of pQCT images 
of trabecular bone in human proximal tibia of different osteoporotic stages. 
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Fig. 7. Determinism DET (A), laminarity LAM (B), mean line length L (C) and 
trapping time TT (D) obtained from four-dimensional RPs constructed from pQCT 
images of trabecular bone in human proximal tibiae with different degree of osteo- 
porosis. 
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Fig. 8. Recurrence rate RR (A) and laminarity LAM (B) compared with trabecular 
network index TNI and structure complexity index (SCI). 



6. 1 Conclusions 



A generalisation of the method of recurrence plots (RPs) and recurrence quan- 
tification analysis (RQA) for the application to higher-dimensional spatial data 
has been proposed here. This new method can be used for 2D image analysis, 
in particular to reveal and quantify recurrent structures in 2D images. Apply- 
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Table 3 

Rank correlation coefficients R for recurrence quantification measures, BMD and 
structural measures of complexity (only significant values are shown). 

2D-RQA BMD 5„ SCI TNI 



RR -0.94 -0.92 -0.91 0.84 

DET -0.65 -0.58 -0.61 0.61 

LAM -0.78 -0.73 -0.75 0.72 

L _ _ _ _ 

TT -0.57 -0.51 - 0.49 



ing this method on model images, we have shown that it is able to distinguish 
typical spatial structures by means of recurrences. As a first application, we 
have used this method for the comparison of CT images of human proximal 
tibia with different degree of osteoporosis. We have found a clear relationship 
between some of the proposed RQA measures and the complexity and ho- 
mogeneity of the trabecular structure. Moreover, this approach can be easily 
extended to higher dimensions, e.g., for 3D analysis of micro-CT images of 
human bone. This approach will be the base for the further development of 
methods for the assessment of structural alteration in trabecular bone with 
osteoporosis in patients on Earth or in space flying personnel in microgravity 
conditions. 
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